BST213: Final exam

Install needed packages

# install.packages("gdata")
library(gdata)
## gdata: Unable to locate valid perl interpreter
## gdata: 
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata: 
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.
## 
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.
## 
## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.
## 
## Attaching package: 'gdata'
## The following object is masked from 'package:stats':
## 
##     nobs
## The following object is masked from 'package:utils':
## 
##     object.size
## The following object is masked from 'package:base':
## 
##     startsWith
# install.packages("car")
library(car)
## Warning: package 'car' was built under R version 3.5.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.5.2
# install.packages("lubridate")
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
# install.packages("gmodels")
library(gmodels)

# install.packages("RQuantLib")
library(RQuantLib)
## Warning: package 'RQuantLib' was built under R version 3.5.3
# install.packages("e1071")
library(e1071)
## Warning: package 'e1071' was built under R version 3.5.3
# install.packages("ggplot2")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.3
# intall.packages("plotly")
library(plotly)
## Warning: package 'plotly' was built under R version 3.5.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# install.packages("survival")
library(survival)

Load database

# Load database
pSERG <- read.csv("D:\\BST213 APPLIED REGRESSION FOR CLINICAL RESEARCH\\FINAL EXAM\\pSERG.csv")

Data cleaning and data transformation

# Transform date of status epilepticus into date format
pSERG$DATESEIZURE <- as.Date(pSERG$DATESEIZURE, format = "%m/%d/%Y")

# Order by date of status epilepticus
pSERG <- pSERG[order(pSERG$PATIENT_LABEL, pSERG$DATESEIZURE), ]

# Delete duplicate episodes from the same patient
pSERG <- pSERG[!duplicated(pSERG$PATIENT_LABEL), ]

# Delete patients with unknown age
pSERG$G_T_STTS_PLPTCS_EPISODE_MONTHS <- as.numeric(as.character(pSERG$G_T_STTS_PLPTCS_EPISODE_MONTHS))
## Warning: NAs introduced by coercion
pSERG$G_T_STTS_PLPTCUS_EPISODE_YEARS <- as.numeric(as.character(pSERG$G_T_STTS_PLPTCUS_EPISODE_YEARS))
## Warning: NAs introduced by coercion
pSERG$ageyears <- pSERG$G_T_STTS_PLPTCUS_EPISODE_YEARS + (pSERG$G_T_STTS_PLPTCS_EPISODE_MONTHS/12)
pSERG <- pSERG[complete.cases(pSERG[ ,"ageyears"]), ]

# Delete patients with unknown sex
pSERG <- pSERG[which(pSERG$SEX == "male" | pSERG$SEX == "female"), ]
pSERG$SEX <- droplevels(pSERG$SEX)

# Delete patients with control (non-refractory) SE 
pSERG <- pSERG[which(pSERG$SE_GROUP == "refractory_case"), ]
pSERG$SE_GROUP <- droplevels(pSERG$SE_GROUP)

# Delete patients with unknown hospital onset
pSERG <- pSERG[which(pSERG$HOSPITALONSET == "no" | pSERG$HOSPITALONSET == "yes"), ]
pSERG$HOSPITALONSET <- droplevels(pSERG$HOSPITALONSET)

# Transform BZDTIME.0 to numeric
pSERG$BZDTIME.0 <- as.numeric(as.character(pSERG$BZDTIME.0))
## Warning: NAs introduced by coercion
# Delete patients with unknown time to first BZD
pSERG <- pSERG[complete.cases(pSERG[ , "BZDTIME.0"]), ]

# Transform AEDTIME.0 to numeric
pSERG$AEDTIME.0 <- as.numeric(as.character(pSERG$AEDTIME.0))
## Warning: NAs introduced by coercion
# Delete patients with unknown time to first non-BZD-AED
pSERG <- pSERG[complete.cases(pSERG[ , "AEDTIME.0"]), ]

# Delete patients with unknown type of SE (continuous vs intermittent)
pSERG <- pSERG[which(pSERG$TYPESTATUS == "continuous" | pSERG$TYPESTATUS == "intermittent"), ]
pSERG$TYPESTATUS <- droplevels(pSERG$TYPESTATUS)

# Create convulsive duration in minutes and eliminate patients with unknown convulsive duration
pSERG$CONVULSIVEDURATION <- as.numeric(as.character(pSERG$CONVULSIVEDURATION))
## Warning: NAs introduced by coercion
pSERG$CONVULSIVEmin <- pSERG$CONVULSIVEDURATION
pSERG$CONVULSIVEhr <- pSERG$CONVULSIVEDURATION * 60
pSERG$convulsivedurationmin <- ifelse(pSERG$CONVULSIVEDURATIONUNITS=="min", pSERG$CONVULSIVEmin, pSERG$CONVULSIVEhr)
pSERG <- pSERG[complete.cases(pSERG[ , "convulsivedurationmin"]), ]

# Delete patients with unknown race
pSERG <- pSERG[complete.cases(pSERG[ , "RACE"]), ]
pSERG$RACE <- droplevels(pSERG$RACE)

# Delete patients with unknown or nonsensical time of SE
pSERG$TIMESEIZURE_HOURS <- as.numeric(as.character(pSERG$TIMESEIZURE_HOURS))
## Warning: NAs introduced by coercion
pSERG <- pSERG[complete.cases(pSERG[ ,"TIMESEIZURE_HOURS"]), ]
pSERG$TIMESEIZURE_HOURS <- ifelse(pSERG$TIMESEIZURE_HOURS >= 100, pSERG$TIMESEIZURE_HOURS/100, pSERG$TIMESEIZURE_HOURS)
# Create variable day/night
pSERG$day <- ifelse(pSERG$TIMESEIZURE_HOURS >= 8 & pSERG$TIMESEIZURE_HOURS < 20, 1, 0)
# Delete patients with time of SE onset
pSERG <- pSERG[complete.cases(pSERG[ , "day"]), ]


###############VARIABLE CREATION#########################
# Divide race into White and non-white
pSERG$white <- ifelse(pSERG$RACE == "white", 1, 0)

# Create variable delay
pSERG$delay[grepl("delay", pSERG$PAST)] <- 1
pSERG$delay[!grepl("delay", pSERG$PAST)] <- 0

# Create variable cerebral palsy
pSERG$palsy[grepl("palsy", pSERG$PAST)] <- 1
pSERG$palsy[!grepl("palsy", pSERG$PAST)] <- 0

# Create variable cerebral palsy or delay
pSERG$delay_or_cerebralpalsy <- ifelse(pSERG$delay==1 | pSERG$palsy==1, 1, 0)

# Create variable none (no neurological comorbidities)
pSERG$none[grepl("none", pSERG$PAST)] <- 1
pSERG$none[!grepl("none", pSERG$PAST)] <- 0

# Create variable prior epilepsy
pSERG$priorepilepsy[grepl("epi",pSERG$PAST)] <- 1
pSERG$priorepilepsy[!grepl("epi",pSERG$PAST)] <- 0

# Create variable prior status
pSERG$status[grepl("status",pSERG$PAST)] <- 1
pSERG$status[!grepl("status",pSERG$PAST)] <- 0

# Create variable of at least one continuous infusion
pSERG$CI <- ifelse(!(pSERG$CONTMED.0==""), 1, 0)

# Transform CI time into numeric
pSERG$CONTTIME.0 <- as.numeric(as.character(pSERG$CONTTIME.0))
## Warning: NAs introduced by coercion
# Create ICU stay in days
pSERG$ICU_DURATION <- as.numeric(as.character(pSERG$ICU_DURATION))
## Warning: NAs introduced by coercion
pSERG$ICUhours <- pSERG$ICU_DURATION/24
pSERG$ICUdays <- pSERG$ICU_DURATION
pSERG$ICUdurationdays <- ifelse(pSERG$ICU_UNITS=="days", pSERG$ICUdays, pSERG$ICUhours)

# Transform EMS arrival to numeric
pSERG$EMSARRIVAL <- as.numeric(as.character(pSERG$EMSARRIVAL))
## Warning: NAs introduced by coercion
# Transform time to hospital arrival to numeric
pSERG$HOSPITALARRIVAL <- as.numeric(as.character(pSERG$HOSPITALARRIVAL))
## Warning: NAs introduced by coercion
# Create variable BZD before hospital arrival
pSERG$BZDbeforehospital <- ifelse(pSERG$BZDTIME.0 < pSERG$HOSPITALARRIVAL, 1, 0)

# Reclasify etiology
pSERG$etiology2 <- recode(pSERG$ETIOLOGY, 
                          "'genetic' = 'genetic'; 
                          'metabolic'= 'metabolic'; 
                          'other' = 'other'; 
                          'structural' = 'structural'; 
                          'unknown' = 'unknown'; 
                          '' = 'unknown'")

# Structural etiology
pSERG$structuraletiology <- ifelse(pSERG$etiology2 == "structural", 1, 0)

# Create variable early in the academic year
pSERG$dateSE <- as.POSIXct(pSERG$DATESEIZURE, format = "%m/%d/%Y")
pSERG$monthSE <- month(pSERG$dateSE, label = FALSE)
pSERG$earlyacademicyear <- ifelse(pSERG$monthSE >= 7 & pSERG$monthSE <= 12, 1, 0)

# Create variable awareness of delays in time to treatment
pSERG$yearSE <- year(pSERG$dateSE)
pSERG$awareness <- ifelse(pSERG$yearSE >= 2015, 1, 0)

# Number of days between the beginning and end of the study
days_of_study <- seq(as.Date("2011-6-1"), as.Date("2019-5-1"), by="days")

# Proportion of holidays in all the days of the study
holiday_days_of_study <- as.factor(ifelse(isBusinessDay(calendar = "UnitedStates", dates = days_of_study), 0, 1))

# Create variable weekends and holidays
pSERG$holiday <- as.factor(ifelse(isBusinessDay(calendar = "UnitedStates", dates = pSERG$DATESEIZURE), 0, 1))

Demographics

# Proportion of patients with ID
CrossTable(pSERG$delay_or_cerebralpalsy)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       149 |       151 | 
##           |     0.497 |     0.503 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Age
nobs(pSERG$ageyears)
## [1] 300
summary(pSERG$ageyears)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.08333  1.25000  3.83333  5.70721  9.25000 20.74167
sd(pSERG$ageyears)
## [1] 5.171068
# Age in patients with ID
nobs(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ageyears)
## [1] 151
summary(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ageyears)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.08333  2.20833  4.91667  6.37369  9.50000 20.74167
sd(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ageyears)
## [1] 5.069913
# Age in patients without ID
nobs(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ageyears)
## [1] 149
summary(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ageyears)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.125   1.000   2.500   5.032   8.917  19.306
sd(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ageyears)
## [1] 5.201736
# Sex
CrossTable(pSERG$SEX)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |    female |      male | 
##           |-----------|-----------|
##           |       128 |       172 | 
##           |     0.427 |     0.573 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Sex in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$SEX)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |    female |      male | 
##           |-----------|-----------|
##           |        68 |        83 | 
##           |     0.450 |     0.550 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Sex in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$SEX)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |    female |      male | 
##           |-----------|-----------|
##           |        60 |        89 | 
##           |     0.403 |     0.597 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Race
CrossTable(pSERG$white)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       110 |       190 | 
##           |     0.367 |     0.633 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Race in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$white)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        52 |        99 | 
##           |     0.344 |     0.656 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Race in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$white)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        58 |        91 | 
##           |     0.389 |     0.611 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Prior epilepsy
CrossTable(pSERG$priorepilepsy)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       156 |       144 | 
##           |     0.520 |     0.480 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Prior epilepsy in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$priorepilepsy)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        40 |       111 | 
##           |     0.265 |     0.735 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Prior epilepsy in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$priorepilepsy)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       116 |        33 | 
##           |     0.779 |     0.221 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Structural etiology
CrossTable(pSERG$structuraletiology)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       222 |        78 | 
##           |     0.740 |     0.260 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Etiology in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$structuraletiology)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       111 |        40 | 
##           |     0.735 |     0.265 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Etiology in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$structuraletiology)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       111 |        38 | 
##           |     0.745 |     0.255 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Hospital onset
CrossTable(pSERG$HOSPITALONSET)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |        no |       yes | 
##           |-----------|-----------|
##           |       207 |        93 | 
##           |     0.690 |     0.310 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Hospital onset in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$HOSPITALONSET)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |        no |       yes | 
##           |-----------|-----------|
##           |       107 |        44 | 
##           |     0.709 |     0.291 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Hospital onset in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$HOSPITALONSET)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |        no |       yes | 
##           |-----------|-----------|
##           |       100 |        49 | 
##           |     0.671 |     0.329 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Type of status epilepticus
CrossTable(pSERG$TYPESTATUS)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##              |   continuous | intermittent | 
##              |--------------|--------------|
##              |           97 |          203 | 
##              |        0.323 |        0.677 | 
##              |--------------|--------------|
## 
## 
## 
## 
# Type of SE in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$TYPESTATUS)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##              |   continuous | intermittent | 
##              |--------------|--------------|
##              |           52 |           99 | 
##              |        0.344 |        0.656 | 
##              |--------------|--------------|
## 
## 
## 
## 
# Type of SE in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$TYPESTATUS)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##              |   continuous | intermittent | 
##              |--------------|--------------|
##              |           45 |          104 | 
##              |        0.302 |        0.698 | 
##              |--------------|--------------|
## 
## 
## 
## 
# Day versus night
CrossTable(pSERG$day)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       126 |       174 | 
##           |     0.420 |     0.580 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Day versus night in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$day)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        71 |        80 | 
##           |     0.470 |     0.530 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Day versus night in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$day)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        55 |        94 | 
##           |     0.369 |     0.631 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Weekday versus weekend/holiday
CrossTable(pSERG$holiday)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       208 |        92 | 
##           |     0.693 |     0.307 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Weekday versus weekend/holiday in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$holiday)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       109 |        42 | 
##           |     0.722 |     0.278 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Weekday versus weekend/holiday in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$holiday)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        99 |        50 | 
##           |     0.664 |     0.336 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Season
CrossTable(pSERG$earlyacademicyear)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       160 |       140 | 
##           |     0.533 |     0.467 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Season in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$earlyacademicyear)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        77 |        74 | 
##           |     0.510 |     0.490 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Season in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$earlyacademicyear)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |        83 |        66 | 
##           |     0.557 |     0.443 | 
##           |-----------|-----------|
## 
## 
## 
## 
# Convulsive duration
nobs(pSERG$convulsivedurationmin)
## [1] 300
summary(pSERG$convulsivedurationmin)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##      0.0     60.0    124.5   2184.0    281.5 172800.0
# Convulsive duration in patients with ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$convulsivedurationmin)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    75.0   140.0   984.6   300.0 37440.0
# Convulsive duration in patients without ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$convulsivedurationmin)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0      53     120    3400     245  172800
# Length of stay
nobs(pSERG$ICUdurationdays)
## [1] 287
summary(pSERG$ICUdurationdays)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00    2.00    4.00   12.08   10.94  180.00      13
# Length of stay in patients with ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ICUdurationdays)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   2.000   4.000   8.532   9.750  73.000       6
# Length of stay in patients without ID
summary(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ICUdurationdays)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   2.000   4.083  15.695  12.000 180.000       7
# Mortality
CrossTable(pSERG$ALIVE)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |           |        no |       yes | 
##           |-----------|-----------|-----------|
##           |         2 |        10 |       288 | 
##           |     0.007 |     0.033 |     0.960 | 
##           |-----------|-----------|-----------|
## 
## 
## 
## 
# Febrile in patients in patients with ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==1, ]$ALIVE)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  151 
## 
##  
##           |           |        no |       yes | 
##           |-----------|-----------|-----------|
##           |         1 |         3 |       147 | 
##           |     0.007 |     0.020 |     0.974 | 
##           |-----------|-----------|-----------|
## 
## 
## 
## 
# Febrile in patients without ID
CrossTable(pSERG[pSERG$delay_or_cerebralpalsy==0, ]$ALIVE)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  149 
## 
##  
##           |           |        no |       yes | 
##           |-----------|-----------|-----------|
##           |         1 |         7 |       141 | 
##           |     0.007 |     0.047 |     0.946 | 
##           |-----------|-----------|-----------|
## 
## 
## 
## 

Univariate analysis time to first BZD

## Outcome: time to first treatment
summary(pSERG$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   17.00   65.53   45.00 1440.00
sd(pSERG$BZDTIME.0)
## [1] 162.3031
skewness(pSERG$BZDTIME.0)
## [1] 5.255733
kurtosis(pSERG$BZDTIME.0)
## [1] 33.28925
plot_firstBZD <- ggplot(aes(x=BZDTIME.0), data=pSERG) + geom_density(color="green", fill="green") +  geom_vline(xintercept=mean(pSERG$BZDTIME.0), color="red", lwd=1) +  geom_vline(xintercept=median(pSERG$BZDTIME.0), color="blue", lwd=1) + theme_classic()
plot_firstBZD

ggplotly(plot_firstBZD)
# Time to first treatment
survfit(Surv(pSERG$BZDTIME.0) ~ 1)
## Call: survfit(formula = Surv(pSERG$BZDTIME.0) ~ 1)
## 
##       n  events  median 0.95LCL 0.95UCL 
##     300     300      17      14      20
plot(survfit(Surv(pSERG$BZDTIME.0) ~ 1), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("purple4"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")

# Time to first BZD depending on ID
summary(pSERG[which(pSERG$delay_or_cerebralpalsy == 1), ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   15.00   61.91   48.50 1264.00
summary(pSERG[which(pSERG$delay_or_cerebralpalsy == 0), ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     5.0    20.0    69.2    42.0  1440.0
survdiff(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy, rho=1)
## Call:
## survdiff(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy, 
##     rho = 1)
## 
##                                  N Observed Expected (O-E)^2/E (O-E)^2/V
## pSERG$delay_or_cerebralpalsy=0 149     76.8     78.1    0.0208    0.0649
## pSERG$delay_or_cerebralpalsy=1 151     79.7     78.4    0.0207    0.0649
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.8
pchisq(survdiff(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy, rho=1)$chisq, length(survdiff(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy, rho=1)$n)-1, lower.tail = FALSE)
## [1] 0.7988799
# Figure time to first BZD by ID
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("ID", "No ID"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by age
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$ageyears))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$ageyears)
## 
##   n= 300, number of events= 300 
## 
##                     coef exp(coef)  se(coef)     z Pr(>|z|)
## pSERG$ageyears 0.0004671 1.0004672 0.0108028 0.043    0.966
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## pSERG$ageyears         1     0.9995    0.9795     1.022
## 
## Concordance= 0.518  (se = 0.021 )
## Rsquare= 0   (max possible= 1 )
## Likelihood ratio test= 0  on 1 df,   p=1
## Wald test            = 0  on 1 df,   p=1
## Score (logrank) test = 0  on 1 df,   p=1
# Time to first treatment by sex
summary(pSERG[pSERG$SEX=="male", ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   20.00   54.38   45.00  720.00
sd(pSERG[pSERG$SEX=="male", ]$BZDTIME.0)
## [1] 107.4961
summary(pSERG[pSERG$SEX=="female", ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   15.00   80.52   46.00 1440.00
sd(pSERG[pSERG$SEX=="female", ]$BZDTIME.0)
## [1] 214.6275
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$SEX))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$SEX)
## 
##   n= 300, number of events= 300 
## 
##                  coef exp(coef) se(coef)    z Pr(>|z|)
## pSERG$SEXmale 0.04609   1.04717  0.11812 0.39    0.696
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## pSERG$SEXmale     1.047      0.955    0.8308      1.32
## 
## Concordance= 0.491  (se = 0.018 )
## Rsquare= 0.001   (max possible= 1 )
## Likelihood ratio test= 0.15  on 1 df,   p=0.7
## Wald test            = 0.15  on 1 df,   p=0.7
## Score (logrank) test = 0.15  on 1 df,   p=0.7
# Figure time to first BZD by sex
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$SEX), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=levels(pSERG$SEX),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by race
summary(pSERG[pSERG$white==1, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   19.00   67.62   47.75 1440.00
sd(pSERG[pSERG$white==1, ]$BZDTIME.0)
## [1] 157.3164
summary(pSERG[pSERG$white==0, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   15.00   61.93   39.50 1264.00
sd(pSERG[pSERG$white==0, ]$BZDTIME.0)
## [1] 171.2509
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$white))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$white)
## 
##   n= 300, number of events= 300 
## 
##                coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$white -0.1113    0.8947   0.1203 -0.925    0.355
## 
##             exp(coef) exp(-coef) lower .95 upper .95
## pSERG$white    0.8947      1.118    0.7068     1.133
## 
## Concordance= 0.516  (se = 0.017 )
## Rsquare= 0.003   (max possible= 1 )
## Likelihood ratio test= 0.85  on 1 df,   p=0.4
## Wald test            = 0.86  on 1 df,   p=0.4
## Score (logrank) test = 0.86  on 1 df,   p=0.4
# Figure time to first BZD by race
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$white), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("non-white", "white"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by prior epilepsy
summary(pSERG[pSERG$priorepilepsy==1, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   15.00   59.78   51.25  660.00
sd(pSERG[pSERG$priorepilepsy==1, ]$BZDTIME.0)
## [1] 119.3409
summary(pSERG[pSERG$priorepilepsy==0, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    7.00   20.00   70.85   40.50 1440.00
sd(pSERG[pSERG$priorepilepsy==0, ]$BZDTIME.0)
## [1] 193.9492
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy)
## 
##   n= 300, number of events= 300 
## 
##                       coef exp(coef) se(coef)     z Pr(>|z|)
## pSERG$priorepilepsy 0.0274    1.0278   0.1166 0.235    0.814
## 
##                     exp(coef) exp(-coef) lower .95 upper .95
## pSERG$priorepilepsy     1.028      0.973    0.8177     1.292
## 
## Concordance= 0.511  (se = 0.018 )
## Rsquare= 0   (max possible= 1 )
## Likelihood ratio test= 0.06  on 1 df,   p=0.8
## Wald test            = 0.06  on 1 df,   p=0.8
## Score (logrank) test = 0.06  on 1 df,   p=0.8
# Figure time to first BZD by prior epilepsy
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("no prior epilepsy", "prior epilepsy"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by structural etiology
summary(pSERG[pSERG$structuraletiology==1, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   15.50   47.74   45.00  460.00
sd(pSERG[pSERG$structuraletiology==1, ]$BZDTIME.0)
## [1] 81.15584
summary(pSERG[pSERG$structuraletiology==0, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   17.00   71.78   45.00 1440.00
sd(pSERG[pSERG$structuraletiology==0, ]$BZDTIME.0)
## [1] 182.1918
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$structuraletiology))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$structuraletiology)
## 
##   n= 300, number of events= 300 
## 
##                             coef exp(coef) se(coef)     z Pr(>|z|)
## pSERG$structuraletiology 0.09882   1.10387  0.13280 0.744    0.457
## 
##                          exp(coef) exp(-coef) lower .95 upper .95
## pSERG$structuraletiology     1.104     0.9059    0.8509     1.432
## 
## Concordance= 0.509  (se = 0.016 )
## Rsquare= 0.002   (max possible= 1 )
## Likelihood ratio test= 0.55  on 1 df,   p=0.5
## Wald test            = 0.55  on 1 df,   p=0.5
## Score (logrank) test = 0.55  on 1 df,   p=0.5
# Figure time to first BZD by structural etiology
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$structuraletiology), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("no structural etiology", "structural etiology"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by hospital onset
summary(pSERG[pSERG$HOSPITALONSET=="no", ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    7.00   20.00   69.13   55.00 1264.00
sd(pSERG[pSERG$HOSPITALONSET=="no", ]$BZDTIME.0)
## [1] 155.9042
summary(pSERG[pSERG$HOSPITALONSET=="yes", ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00    9.00   57.53   25.00 1440.00
sd(pSERG[pSERG$HOSPITALONSET=="yes", ]$BZDTIME.0)
## [1] 176.3347
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET)
## 
##   n= 300, number of events= 300 
## 
##                          coef exp(coef) se(coef)     z Pr(>|z|)   
## pSERG$HOSPITALONSETyes 0.3299    1.3908   0.1260 2.619  0.00882 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                        exp(coef) exp(-coef) lower .95 upper .95
## pSERG$HOSPITALONSETyes     1.391      0.719     1.087      1.78
## 
## Concordance= 0.563  (se = 0.016 )
## Rsquare= 0.022   (max possible= 1 )
## Likelihood ratio test= 6.56  on 1 df,   p=0.01
## Wald test            = 6.86  on 1 df,   p=0.009
## Score (logrank) test = 6.92  on 1 df,   p=0.009
# Figure time to first BZD by hospital onset
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("out of the hospital", "in the hospital"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by type of status epilepticus
summary(pSERG[pSERG$TYPESTATUS=="intermittent", ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   20.00   85.71   60.00 1440.00
sd(pSERG[pSERG$TYPESTATUS=="intermittent", ]$BZDTIME.0)
## [1] 193.1872
summary(pSERG[pSERG$TYPESTATUS=="continuous", ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    6.00   14.00   23.31   28.00  180.00
sd(pSERG[pSERG$TYPESTATUS=="continuous", ]$BZDTIME.0)
## [1] 29.20201
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS)
## 
##   n= 300, number of events= 300 
## 
##                                 coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$TYPESTATUSintermittent -0.4991    0.6071   0.1296 -3.852 0.000117
##                                 
## pSERG$TYPESTATUSintermittent ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$TYPESTATUSintermittent    0.6071      1.647    0.4709    0.7826
## 
## Concordance= 0.536  (se = 0.017 )
## Rsquare= 0.046   (max possible= 1 )
## Likelihood ratio test= 14.12  on 1 df,   p=2e-04
## Wald test            = 14.84  on 1 df,   p=1e-04
## Score (logrank) test = 15.13  on 1 df,   p=1e-04
# Figure time to first BZD by type of status epilepticus
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("continuous", "intermittent"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by day versus night
summary(pSERG[pSERG$day==1, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   16.00   54.95   41.75 1132.00
sd(pSERG[pSERG$day==1, ]$BZDTIME.0)
## [1] 126.6443
summary(pSERG[pSERG$day==0, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   20.00   80.14   50.00 1440.00
sd(pSERG[pSERG$day==0, ]$BZDTIME.0)
## [1] 201.1023
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$day))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$day)
## 
##   n= 300, number of events= 300 
## 
##             coef exp(coef) se(coef)     z Pr(>|z|)
## pSERG$day 0.1230    1.1308   0.1182 1.041    0.298
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## pSERG$day     1.131     0.8843    0.8971     1.426
## 
## Concordance= 0.509  (se = 0.018 )
## Rsquare= 0.004   (max possible= 1 )
## Likelihood ratio test= 1.09  on 1 df,   p=0.3
## Wald test            = 1.08  on 1 df,   p=0.3
## Score (logrank) test = 1.08  on 1 df,   p=0.3
# Figure time to first BZD by day versus night
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$day), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("night", "day"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by weekday versus weekend/holiday
summary(pSERG[pSERG$holiday==0, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    6.00   20.00   73.67   50.00 1264.00
sd(pSERG[pSERG$holiday==0, ]$BZDTIME.0)
## [1] 163.9659
summary(pSERG[pSERG$holiday==1, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   10.50   47.14   31.25 1440.00
sd(pSERG[pSERG$holiday==1, ]$BZDTIME.0)
## [1] 157.8114
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$holiday))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$holiday)
## 
##   n= 300, number of events= 300 
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)  
## pSERG$holiday1 0.2851    1.3299   0.1265 2.254   0.0242 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## pSERG$holiday1      1.33     0.7519     1.038     1.704
## 
## Concordance= 0.541  (se = 0.016 )
## Rsquare= 0.016   (max possible= 1 )
## Likelihood ratio test= 4.89  on 1 df,   p=0.03
## Wald test            = 5.08  on 1 df,   p=0.02
## Score (logrank) test = 5.11  on 1 df,   p=0.02
# Figure time to first BZD by weekday versus weekend/holiday
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$holiday), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("weekday", "weekend/holiday"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

# Time to first treatment by season in the academic year
summary(pSERG[pSERG$earlyacademicyear==1, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   15.00   55.06   35.75 1264.00
sd(pSERG[pSERG$earlyacademicyear==1, ]$BZDTIME.0)
## [1] 140.9172
summary(pSERG[pSERG$earlyacademicyear==0, ]$BZDTIME.0)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    6.75   20.00   74.70   55.00 1440.00
sd(pSERG[pSERG$earlyacademicyear==0, ]$BZDTIME.0)
## [1] 178.8734
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$earlyacademicyear))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$earlyacademicyear)
## 
##   n= 300, number of events= 300 
## 
##                           coef exp(coef) se(coef)     z Pr(>|z|)
## pSERG$earlyacademicyear 0.1752    1.1915   0.1162 1.508    0.131
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## pSERG$earlyacademicyear     1.192     0.8393    0.9489     1.496
## 
## Concordance= 0.529  (se = 0.018 )
## Rsquare= 0.008   (max possible= 1 )
## Likelihood ratio test= 2.26  on 1 df,   p=0.1
## Wald test            = 2.28  on 1 df,   p=0.1
## Score (logrank) test = 2.28  on 1 df,   p=0.1
# Figure time to first BZD by season in the academic year
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$earlyacademicyear), fun = "event", 
     conf.int = FALSE, xlim = c(0,60), col = c("violetred3", "seagreen3"), lwd = 3, 
     cex.axis = 1.3, cex.lab = 1.5,
     frame.plot=FALSE,
     xlab= "Time to first BZD (min)", ylab= "Cum. prob. having received first BZD")
legend("topleft", legend=c("January-June", "July-December"),
  col=c("violetred3", "seagreen3"),
  lty=1,
  lwd=3,
  horiz=FALSE,
  bty='n',
  cex=1.3)

Multivariable model time to first BZD

# Initial model with the candidate variables with a p-value less than 0.25 on multivariable analysis:
#hospital onset, type of SE, holiday, and season in the academic year
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$earlyacademicyear))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$earlyacademicyear)
## 
##   n= 300, number of events= 300 
## 
##                                 coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.0667    1.0690   0.1172  0.569  0.56938
## pSERG$HOSPITALONSETyes        0.3552    1.4265   0.1283  2.769  0.00562
## pSERG$TYPESTATUSintermittent -0.5736    0.5635   0.1320 -4.345  1.4e-05
## pSERG$holiday1                0.3075    1.3600   0.1296  2.372  0.01769
## pSERG$earlyacademicyear       0.1300    1.1388   0.1173  1.108  0.26793
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$earlyacademicyear         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0690     0.9355    0.8495    1.3451
## pSERG$HOSPITALONSETyes          1.4265     0.7010    1.1094    1.8342
## pSERG$TYPESTATUSintermittent    0.5635     1.7746    0.4350    0.7299
## pSERG$holiday1                  1.3600     0.7353    1.0549    1.7535
## pSERG$earlyacademicyear         1.1388     0.8781    0.9049    1.4332
## 
## Concordance= 0.601  (se = 0.021 )
## Rsquare= 0.096   (max possible= 1 )
## Likelihood ratio test= 30.28  on 5 df,   p=1e-05
## Wald test            = 31.37  on 5 df,   p=8e-06
## Score (logrank) test = 31.82  on 5 df,   p=6e-06
# Model without season in the academic year because p-value >0.1 in the multivariable model and it is not 
#a confounder because it does not change the coefficient of intellectual disability by more than 20%
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday)
## 
##   n= 300, number of events= 300 
## 
##                                  coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.07959   1.08285  0.11669  0.682  0.49516
## pSERG$HOSPITALONSETyes        0.35566   1.42712  0.12855  2.767  0.00566
## pSERG$TYPESTATUSintermittent -0.58428   0.55751  0.13198 -4.427 9.56e-06
## pSERG$holiday1                0.31743   1.37360  0.12954  2.450  0.01427
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0828     0.9235    0.8615    1.3611
## pSERG$HOSPITALONSETyes          1.4271     0.7007    1.1093    1.8360
## pSERG$TYPESTATUSintermittent    0.5575     1.7937    0.4304    0.7221
## pSERG$holiday1                  1.3736     0.7280    1.0656    1.7706
## 
## Concordance= 0.599  (se = 0.021 )
## Rsquare= 0.092   (max possible= 1 )
## Likelihood ratio test= 29.06  on 4 df,   p=8e-06
## Wald test            = 30.03  on 4 df,   p=5e-06
## Score (logrank) test = 30.5  on 4 df,   p=4e-06
# Add non-candidate variables to the multivariable model: age (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$ageyears))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$ageyears)
## 
##   n= 300, number of events= 300 
## 
##                                   coef exp(coef)  se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.070733  1.073295  0.119648  0.591  0.55440
## pSERG$HOSPITALONSETyes        0.352977  1.423298  0.128753  2.741  0.00612
## pSERG$TYPESTATUSintermittent -0.591457  0.553520  0.133788 -4.421 9.83e-06
## pSERG$holiday1                0.317058  1.373082  0.129529  2.448  0.01437
## pSERG$ageyears                0.003761  1.003768  0.011263  0.334  0.73842
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$ageyears                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0733     0.9317    0.8489    1.3569
## pSERG$HOSPITALONSETyes          1.4233     0.7026    1.1059    1.8319
## pSERG$TYPESTATUSintermittent    0.5535     1.8066    0.4258    0.7195
## pSERG$holiday1                  1.3731     0.7283    1.0652    1.7699
## pSERG$ageyears                  1.0038     0.9962    0.9819    1.0262
## 
## Concordance= 0.597  (se = 0.021 )
## Rsquare= 0.093   (max possible= 1 )
## Likelihood ratio test= 29.17  on 5 df,   p=2e-05
## Wald test            = 30.09  on 5 df,   p=1e-05
## Score (logrank) test = 30.57  on 5 df,   p=1e-05
# Add non-candidate variables to the multivariable model: sex (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$SEX))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$SEX)
## 
##   n= 300, number of events= 300 
## 
##                                  coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.07884   1.08203  0.11666  0.676  0.49915
## pSERG$HOSPITALONSETyes        0.36019   1.43360  0.12885  2.795  0.00518
## pSERG$TYPESTATUSintermittent -0.57804   0.56100  0.13253 -4.362 1.29e-05
## pSERG$holiday1                0.32484   1.38380  0.13038  2.492  0.01272
## pSERG$SEXmale                 0.05958   1.06139  0.12035  0.495  0.62054
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$SEXmale                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy     1.082     0.9242    0.8609    1.3600
## pSERG$HOSPITALONSETyes           1.434     0.6975    1.1137    1.8455
## pSERG$TYPESTATUSintermittent     0.561     1.7825    0.4327    0.7274
## pSERG$holiday1                   1.384     0.7226    1.0718    1.7867
## pSERG$SEXmale                    1.061     0.9422    0.8384    1.3437
## 
## Concordance= 0.599  (se = 0.021 )
## Rsquare= 0.093   (max possible= 1 )
## Likelihood ratio test= 29.3  on 5 df,   p=2e-05
## Wald test            = 30.3  on 5 df,   p=1e-05
## Score (logrank) test = 30.76  on 5 df,   p=1e-05
# Add non-candidate variables to the multivariable model: white (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$white))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$white)
## 
##   n= 300, number of events= 300 
## 
##                                  coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.08515   1.08888  0.11695  0.728  0.46653
## pSERG$HOSPITALONSETyes        0.36022   1.43364  0.12877  2.797  0.00515
## pSERG$TYPESTATUSintermittent -0.57366   0.56346  0.13272 -4.322 1.54e-05
## pSERG$holiday1                0.31720   1.37328  0.12960  2.448  0.01438
## pSERG$white                  -0.08857   0.91524  0.12149 -0.729  0.46598
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$white                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0889     0.9184    0.8658    1.3694
## pSERG$HOSPITALONSETyes          1.4336     0.6975    1.1139    1.8452
## pSERG$TYPESTATUSintermittent    0.5635     1.7748    0.4344    0.7309
## pSERG$holiday1                  1.3733     0.7282    1.0652    1.7704
## pSERG$white                     0.9152     1.0926    0.7213    1.1613
## 
## Concordance= 0.601  (se = 0.021 )
## Rsquare= 0.094   (max possible= 1 )
## Likelihood ratio test= 29.58  on 5 df,   p=2e-05
## Wald test            = 30.68  on 5 df,   p=1e-05
## Score (logrank) test = 31.13  on 5 df,   p=9e-06
# Add non-candidate variables to the multivariable model: prior epilepsy (enters the model: p-value>0.1 but confounder because it changes the coefficient for ID from 0.07959 to 0.03839, more than 20%)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$priorepilepsy)
## 
##   n= 300, number of events= 300 
## 
##                                  coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.03839   1.03914  0.14013  0.274  0.78411
## pSERG$HOSPITALONSETyes        0.35980   1.43304  0.12878  2.794  0.00521
## pSERG$TYPESTATUSintermittent -0.59333   0.55249  0.13311 -4.457  8.3e-06
## pSERG$holiday1                0.31362   1.36837  0.12964  2.419  0.01556
## pSERG$priorepilepsy           0.07492   1.07780  0.14115  0.531  0.59554
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$priorepilepsy             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0391     0.9623    0.7896    1.3676
## pSERG$HOSPITALONSETyes          1.4330     0.6978    1.1134    1.8445
## pSERG$TYPESTATUSintermittent    0.5525     1.8100    0.4256    0.7172
## pSERG$holiday1                  1.3684     0.7308    1.0613    1.7642
## pSERG$priorepilepsy             1.0778     0.9278    0.8173    1.4213
## 
## Concordance= 0.598  (se = 0.021 )
## Rsquare= 0.093   (max possible= 1 )
## Likelihood ratio test= 29.34  on 5 df,   p=2e-05
## Wald test            = 30.22  on 5 df,   p=1e-05
## Score (logrank) test = 30.72  on 5 df,   p=1e-05
# Add non-candidate variables to the multivariable model: structural etiology (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy + pSERG$structuraletiology))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$priorepilepsy + pSERG$structuraletiology)
## 
##   n= 300, number of events= 300 
## 
##                                  coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.04230   1.04320  0.13998  0.302  0.76253
## pSERG$HOSPITALONSETyes        0.34645   1.41404  0.12987  2.668  0.00764
## pSERG$TYPESTATUSintermittent -0.59226   0.55308  0.13285 -4.458 8.27e-06
## pSERG$holiday1                0.33052   1.39170  0.13144  2.515  0.01192
## pSERG$priorepilepsy           0.07784   1.08095  0.14087  0.553  0.58059
## pSERG$structuraletiology      0.11021   1.11652  0.13568  0.812  0.41660
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$priorepilepsy             
## pSERG$structuraletiology        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0432     0.9586    0.7929    1.3725
## pSERG$HOSPITALONSETyes          1.4140     0.7072    1.0963    1.8239
## pSERG$TYPESTATUSintermittent    0.5531     1.8081    0.4263    0.7176
## pSERG$holiday1                  1.3917     0.7185    1.0756    1.8006
## pSERG$priorepilepsy             1.0809     0.9251    0.8201    1.4247
## pSERG$structuraletiology        1.1165     0.8956    0.8558    1.4566
## 
## Concordance= 0.598  (se = 0.021 )
## Rsquare= 0.095   (max possible= 1 )
## Likelihood ratio test= 29.99  on 6 df,   p=4e-05
## Wald test            = 30.97  on 6 df,   p=3e-05
## Score (logrank) test = 31.42  on 6 df,   p=2e-05
# Add non-candidate variables to the multivariable model: day versus night (does not enter the model: p-value>0.1 and non-confounder)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy + pSERG$day))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$priorepilepsy + pSERG$day)
## 
##   n= 300, number of events= 300 
## 
##                                  coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.03354   1.03411  0.14046  0.239  0.81129
## pSERG$HOSPITALONSETyes        0.36118   1.43502  0.12846  2.812  0.00493
## pSERG$TYPESTATUSintermittent -0.58839   0.55522  0.13319 -4.418 9.98e-06
## pSERG$holiday1                0.31087   1.36461  0.12954  2.400  0.01641
## pSERG$priorepilepsy           0.09301   1.09747  0.14262  0.652  0.51431
## pSERG$day                     0.11949   1.12692  0.11960  0.999  0.31774
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$priorepilepsy             
## pSERG$day                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0341     0.9670    0.7852    1.3618
## pSERG$HOSPITALONSETyes          1.4350     0.6969    1.1156    1.8459
## pSERG$TYPESTATUSintermittent    0.5552     1.8011    0.4277    0.7208
## pSERG$holiday1                  1.3646     0.7328    1.0586    1.7590
## pSERG$priorepilepsy             1.0975     0.9112    0.8298    1.4514
## pSERG$day                       1.1269     0.8874    0.8914    1.4246
## 
## Concordance= 0.6  (se = 0.021 )
## Rsquare= 0.096   (max possible= 1 )
## Likelihood ratio test= 30.34  on 6 df,   p=3e-05
## Wald test            = 31.4  on 6 df,   p=2e-05
## Score (logrank) test = 31.88  on 6 df,   p=2e-05
# Final multivariable model
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$priorepilepsy)
## 
##   n= 300, number of events= 300 
## 
##                                  coef exp(coef) se(coef)      z Pr(>|z|)
## pSERG$delay_or_cerebralpalsy  0.03839   1.03914  0.14013  0.274  0.78411
## pSERG$HOSPITALONSETyes        0.35980   1.43304  0.12878  2.794  0.00521
## pSERG$TYPESTATUSintermittent -0.59333   0.55249  0.13311 -4.457  8.3e-06
## pSERG$holiday1                0.31362   1.36837  0.12964  2.419  0.01556
## pSERG$priorepilepsy           0.07492   1.07780  0.14115  0.531  0.59554
##                                 
## pSERG$delay_or_cerebralpalsy    
## pSERG$HOSPITALONSETyes       ** 
## pSERG$TYPESTATUSintermittent ***
## pSERG$holiday1               *  
## pSERG$priorepilepsy             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                              exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy    1.0391     0.9623    0.7896    1.3676
## pSERG$HOSPITALONSETyes          1.4330     0.6978    1.1134    1.8445
## pSERG$TYPESTATUSintermittent    0.5525     1.8100    0.4256    0.7172
## pSERG$holiday1                  1.3684     0.7308    1.0613    1.7642
## pSERG$priorepilepsy             1.0778     0.9278    0.8173    1.4213
## 
## Concordance= 0.598  (se = 0.021 )
## Rsquare= 0.093   (max possible= 1 )
## Likelihood ratio test= 29.34  on 5 df,   p=2e-05
## Wald test            = 30.22  on 5 df,   p=1e-05
## Score (logrank) test = 30.72  on 5 df,   p=1e-05
# Test for proportionality 
time.dep.zphfirstBZD <- cox.zph(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy))
time.dep.zphfirstBZD
##                                  rho  chisq       p
## pSERG$delay_or_cerebralpalsy  0.0399  0.507 0.47628
## pSERG$HOSPITALONSETyes       -0.1786  9.636 0.00191
## pSERG$TYPESTATUSintermittent -0.1269  5.220 0.02233
## pSERG$holiday1               -0.0249  0.189 0.66355
## pSERG$priorepilepsy          -0.0433  0.609 0.43513
## GLOBAL                            NA 19.607 0.00148
# Proportional hazards assumtion for ID
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy), fun = "cloglog",
 conf.int=FALSE, col = c("red", "blue"), lwd = 3,
 cex.axis = 1.5, cex.lab = 1.5,
 xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")

# Plot residual plots for ID
plot(time.dep.zphfirstBZD[1],
 lwd = 2,
 col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)

# Proportional hazards assumtion for hospital onset
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$HOSPITALONSET), fun = "cloglog",
 conf.int=FALSE, col = c("red", "blue"), lwd = 3,
 cex.axis = 1.5, cex.lab = 1.5,
 xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")

# Plot residual plots for hospital onset
plot(time.dep.zphfirstBZD[2],
 lwd = 2,
 col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)

# Proportional hazards assumtion for type of SE
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$TYPESTATUS), fun = "cloglog",
 conf.int=FALSE, col = c("red", "blue"), lwd = 3,
 cex.axis = 1.5, cex.lab = 1.5,
 xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")

# Plot residual plots for type of SE
plot(time.dep.zphfirstBZD[3],
 lwd = 2,
 col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)

# Proportional hazards assumtion for holiday
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$holiday), fun = "cloglog",
 conf.int=FALSE, col = c("red", "blue"), lwd = 3,
 cex.axis = 1.5, cex.lab = 1.5,
 xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")

# Plot residual plots for holiday
plot(time.dep.zphfirstBZD[4],
 lwd = 2,
 col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)

# Proportional hazards assumtion for prior epilepsy
plot(survfit(Surv(pSERG$BZDTIME.0) ~ pSERG$priorepilepsy), fun = "cloglog",
 conf.int=FALSE, col = c("red", "blue"), lwd = 3,
 cex.axis = 1.5, cex.lab = 1.5,
 xlab="log[Time to first BZD (min)]", ylab="log[Cum. prob. having received first BZD]")

# Plot residual plots for holiday
plot(time.dep.zphfirstBZD[5],
 lwd = 2,
 col = c("magenta3", "purple"))
abline(h=0, lty = 3, lwd = 2)

Multivariable model time to first BZD with interaction with hospital onset

# Create the variable HOSPITALONSETNUMERIC
pSERG$HOSPITALONSETNUMERIC <- ifelse(pSERG$HOSPITALONSET=="yes", 1, 0)

# Create the variable interactionIDhospitalonset
pSERG$interactionIDhospitalonset <- pSERG$delay_or_cerebralpalsy * pSERG$HOSPITALONSETNUMERIC

# Final multivariable model with interaction with hospital onset (non-significant)
summary(coxph(Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + pSERG$priorepilepsy + pSERG$interactionIDhospitalonset))
## Call:
## coxph(formula = Surv(pSERG$BZDTIME.0) ~ pSERG$delay_or_cerebralpalsy + 
##     pSERG$HOSPITALONSET + pSERG$TYPESTATUS + pSERG$holiday + 
##     pSERG$priorepilepsy + pSERG$interactionIDhospitalonset)
## 
##   n= 300, number of events= 300 
## 
##                                      coef exp(coef) se(coef)      z
## pSERG$delay_or_cerebralpalsy      0.01144   1.01150  0.15663  0.073
## pSERG$HOSPITALONSETyes            0.31039   1.36396  0.18259  1.700
## pSERG$TYPESTATUSintermittent     -0.58623   0.55642  0.13441 -4.362
## pSERG$holiday1                    0.31114   1.36498  0.12982  2.397
## pSERG$priorepilepsy               0.06861   1.07101  0.14205  0.483
## pSERG$interactionIDhospitalonset  0.09982   1.10497  0.25969  0.384
##                                  Pr(>|z|)    
## pSERG$delay_or_cerebralpalsy       0.9418    
## pSERG$HOSPITALONSETyes             0.0891 .  
## pSERG$TYPESTATUSintermittent     1.29e-05 ***
## pSERG$holiday1                     0.0165 *  
## pSERG$priorepilepsy                0.6291    
## pSERG$interactionIDhospitalonset   0.7007    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                  exp(coef) exp(-coef) lower .95 upper .95
## pSERG$delay_or_cerebralpalsy        1.0115     0.9886    0.7441    1.3750
## pSERG$HOSPITALONSETyes              1.3640     0.7332    0.9536    1.9508
## pSERG$TYPESTATUSintermittent        0.5564     1.7972    0.4276    0.7241
## pSERG$holiday1                      1.3650     0.7326    1.0583    1.7605
## pSERG$priorepilepsy                 1.0710     0.9337    0.8107    1.4149
## pSERG$interactionIDhospitalonset    1.1050     0.9050    0.6642    1.8382
## 
## Concordance= 0.596  (se = 0.021 )
## Rsquare= 0.094   (max possible= 1 )
## Likelihood ratio test= 29.49  on 6 df,   p=5e-05
## Wald test            = 30.51  on 6 df,   p=3e-05
## Score (logrank) test = 31.07  on 6 df,   p=2e-05